Equilibration of a Tonks-Girardeau gas following a trap release 
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We study the non-equilibrium dynamics of a Tonks-Girardeau gas released from a parabolic trap to a circle. 
We present the exact analytic solution of the many body dynamics and prove that, for large times and in a prop- 
erly defined thermodynamic limit, the reduced density matrix of any finite subsystem converges to a generalized 
Gibbs ensemble. The equilibration mechanism is expected to be the same for all one-dimensional systems. 

PACS numbers: 



The non-equilibrium dynamics of isolated many body 
quantum systems is currently in a golden age mainly due to 
the experiments on trapped ultra-cold atomic gases [1-8] in 
which it is possible to measure the unitary non-equilibrium 
evolution without any significant coupling to the environment. 
A key question is whether the system relaxes to a stationary 
state, and if it does, how to characterize from first principles its 
physical properties at late times. It is commonly believed that, 
depending on the integrability of the Hamiltonian governing 
the time evolution, the behavior of local observables either can 
be described by an effective thermal distribution or by a gen- 
eralized Gibbs ensemble (GGE), for non-integrable and inte- 
grable systems respectively (see e.g. [9] for a review). While 
this scenario is coiToborated by many investigations [10-29], 
a few studies [30-35] suggest that the behavior could be more 
complicated and in particular can depend on the initial state. 

In a global quantum quench, the initial condition is the 
ground state of a translationally invariant Hamiltonian which 
differs from the one governing the evolution by an experimen- 
tally tunable parameter such as a magnetic field. A differ- 
ent initial condition can be experimentally achieved [7, 8] by 
considering the non-equilibrium dynamics of a gas released 
from a parabolic trapping potential. It has been shown exper- 
imentally that the spreading of correlations is ballistic for an 
integrable system and diffusive for a non-integrable one [8]. 
However, when the gas expands in full space, for infinite time 
the gas clearly reaches zero density (see e.g. [35^0] for a 
theoretical analysis) and it is rather confusing to distinguish 
thermal and GGE states. To circumvent this, Caux and Konik 
[41] have recently developed a new approach based on inte- 
grability to study the release of the Lieb-Liniger Bose gas [42] 
from a parabolic trap not in free space but on a closed circle 
(as sketched in Fig. 1), so that the gas has finite density. It 
has been numerically shown that the time averaged correla- 
tion functions are described by a GGE, apart from finite size 
effects [41]. A preliminary analysis for non-integrable models 
has also been presented [43]. However, while this approach 
permits effectively to calculate time-averaged quantities for 
relatively large systems (the maximum number of particles is 

= 56 [41]), the study of the time evolution is possible but 
much harder and it is difficult to establish whether (and in 
which sense) an infinite time limit exists. 

In order to overcome these limitations, we present here 
a full analytic solution of this non-equilibrium dynamics in 




FIG. 1 : Left: sketch of the trap release dynamic in a circle. Right: 
Color plot of the numerical calculated density evolution for A'^ = 
10, 100, oo (from left to right) at N/L = 1 /2 and cjTV = 5. 



the limit of strong coupling, i.e. in the celebrated Tonks- 
Girardeau regime [44]. We will show that, in a properly de- 
fined thermodynamic (TD) limit, the reduced density matrix 
of any finite subsystem converges for long times to the GGE 
one. This implies that any measurable local observable will 
converge to the GGE predictions. 

The model and quench protocol. We consider a one- 
dimensional Bose gas with delta pairwise interaction and in 
an external parabolic potential with Hamiltonian 



H = - 



I ^ Q2 N ^ 
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where c > is the coupling constant (we set S = to = 1). 
The translationally invariant Lieb-Liniger model is obtained 
for Lu — Q and on a circle of length L with periodic boundary 
conditions (PBC). While Ref. [41] covers numerically arbi- 
trary c, to make an analytic progress we consider the strong- 
coupling limit of impenetrable bosons c — > oo, corresponding 
also to the low density n = N/L ^ 1 regime for any c [42]. 

For a trap release, the initial state is a Tonks Girardeau gas 
confined by a parabolic potential, i.e. the ground state of Eq. 
(1) for a fixed uj. Following [44], the many body wave func- 
tion for the ground state of N impenetrable bosons is 

*B(a;i, ■ ■ • ,xn) = JJsgn(xj - Xi)'^F{xi, ■ ■ ■ ,a;jv), (2) 

i<3 



2 



where ^ p {xi , • • • , xjv) is the ground-state function of N free 
fermions in the parabolic potential, i.e. the Slater determinant 
detj^j Xji^i) with the eigenstates of the harmonic oscillator 
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and Hj{z) the Hermite polynomials. In this fermionic lan- 
guage, the time evolution governed by the Hamitonian (1) 
with w = is obtained by expanding the one-particle states in 
the free- wave basis, i.e. (k ~ 27nn/L) 
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We now make the only crucial physical assumption: we im- 
pose that the space initially occupied by the trapped gas as a 
whole is within the external box of length L, i.e. before the 
quench the PBC are irrelevant for the gas which only "sees" 
the parabolic trap. This condition is what allows us to talk 
about release of the gas and requires the number of particles 
N to be smaller than the first level of the parabolic potential 
that is affected by PBC. In the TD limit, for large quantum 
numbers, |xAf(2;)P is the semiclassical probability density at 
the corresponding energy that tends to zero for > 1/2 with 
t the classical cloud dimension i = 2^j2Njuj. In simpler 
words, this means that the classical extension of the gas in the 
trap £ must be smaller than the box size L. 

To have a well-defined TD limit, we should consider 
iV, L — ^ oo at fixed density n = N /L and, at the same time, 
a; — > with ujN constant (fixed initial density), as in [41]. 
In terms of these quantities the gas release condition £ < L 
reads V Nuj > 2\/2n and the coefficients Ak,,j can be calcu- 
lated extending the integration in Eq. (4) to ±C)o, obtaining 
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Also the infinite time limit should be handled with care. In- 
deed, in this quench, a stationary behavior is possible because 
of the interference of the particles going around the circle L 
many times (see Fig. 1), i.e. to observe a stationary value we 
must require vt ^ L (the speed of sound is u = m 
our normalization). This is very different from equilibration in 
standard global quenches where the time should be such that 
the boundaries are not reached (see e.g. [17]) in order to avoid 
revival effects. In this problem the revival scale is t,- oc and 
so the infinite time limit in which a stationary behavior can be 
achieved is t/L — > oo provided t/L^ 0. The importance 
of the TD and long time limits to get a stationary behavior is 
akeady evident from the time evolution of the density profile 
in Fig. 1. 

The one-particle problem. In fermion language, the time 
dependent many body state is the Slater determinant of the 
time evolved one-particle initial eigenfunctions (i.e. the so- 
lution of the Schrodinger equation idt^j [x, t) = Hp^j (a;, t) 
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FIG. 2: (a,b,c) Time evolution of the density n{x, t) for different 
X / L and sizes. Dashed red lines indicate the equilibration value N/ L 
at infinite time, (d) Density profile for L — 1600 at different reseated 
times t/L. Symbols are the exact dynamics for finite A'^, while full 
black lines are the TD limit. 



with $j {x, 0) = Xj {x) ™d Hp the single particle free Hamil- 
tonian with PBC). These time evolved wave functions can be 
calculated from Eqs. (4) and (5) obtaining 
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where 
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is the time evolved eigenfunction in infinite space which 
agrees with the result in [36]. The boson-fermion mapping 
remains valid for the time-dependent problem [45]. 

Time evolution of the density profile. We start the analysis 
of the many body problem from the density profile n{x, t) = 
J2j \^j {x, P which shows clearly how a non-zero stationary 
value can be achieved in a trap release experiment. From Eq. 
(7) we have for arbitrary time, N, L, u) 
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which, in the TD limit, because of the strongly oscillating 
phase factor, reduces to the diagonal part p = q: 



n(x, t) 
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and can be rewritten in terms of the TD limit of the particle 
density at initial time 710(2;) = {\/2Nlj — uj'^x^)/'k as 
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In Figs. 1 and 2 we show the numerically calculated time de- 
pendent density for finite but large N which perfectly agrees 
with the above TD prediction for any time. 

The two-point fermionic correlator C{x,y]t) = 
{^^ (x, t)^/{y, t)) is given by 



C{x,y;t) 



E 
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The numerical determination of this correlation function for 
finite N is reported in Fig. 3 showing the approach to the 
infinite time limit [46] 



C{x, y;t ^ 00) — 2n 



Ji [V2ujN{x - y)] 
V2ujN{x ~ y) 



(12) 



with Ji{z) the Bessel function. 

The reduced density matrix and the GGE. For a closed sys- 
tem evolving under Hamiltonian dynamics, the existence of a 
stationary state may seem paradoxical because the whole sys- 
tem is always in a pure state and cannot be described by a 
mixed state at infinite time. This 'paradox' is solved in the 
reduced density matrix formalism: given a space interval A, 
the reduced density matrix is pA{t) = Tr Bp{t) where B is 
the complement of A and p{t) = |*(t))(*(t)| is the density 
matrix of the whole system. With some abuse of language, 
we say that a system becomes stationary if, after the TD limit 
is taken for the whole system, the limit pA.oo — lim PA{t) 

exists for any finite A [17]. Furthermore we say that a system 
is described by a statistical ensemble with density matrix p^ 
if the reduced density matrix pa.e = T^J^bPe equals pA,oo- 

For a gas of free fermions, by means of Wick theorem, any 
observable can be obtained from the two-point correlator The 
construction of pA in terms of C{x, y) in continuous space has 
been detailed in [47] (generalizing the lattice approach [48]). 
As a crucial point, the non-local transformation mapping the 
Tonks-Girardeau gas to free fermions is local within any given 
compact subspace, i.e. the bosonic degrees of freedom within 
A can be written only in terms of fermions in A. This is anal- 
ogous to lattice models such as the Ising chain [16-18]. Thus, 
if for finite a;, y, C{x, y,t 00) is described by a statistical 
ensemble, also pA will be and consequently the expectation 
value of any observable local within A. 
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FIG. 3: Snapshots of the correlation Re[C(x,0;t)] at different 
reseated times t/L and sizes. For t/L — the full line is the initial 
correlation in the TD limit, i.e. C{x, 0) = sin[\/ 2ujNx] / tvx valid 
for X <^ L. The full line for t/L — 2, 4 is the stationary value in Eq. 
(12). As time increases, two symmetric peaks are expelled from the 
central region. The inset in (d) shows the evolution for fixed x — 5 
and L — 1600: after the moving peak has been expelled, the corre- 
lation is damped in time and converges to the GGE. 



Because of integrability, it is natural to expect that Eq. (12) 
should be described by a GGE 



(13) 



with {li} a complete set of local integrals of motion and 
Lagrange multipliers fixed by the conditions (5'o|/i|5'o) — 
Tr[pGG£;-^j]' with |^o) the many body initial state. However, 
for free fermions, one can work with the momentum occupa- 
tion modes hk = c^Ck which are non-local integrals of mo- 
tion, but can be written as linear combinations of local inte- 
grals of motion [28]. In the TD limit, the initial values of fik 
are 



(14) 



and zero if the argument of the square root is negative. In 
the GGE we have Ti[pGGEnk] — {e^'' + and equating 
the two, the are derived. It is now straightforward to show 
that C{x, y) in the GGE equals the infinite time limit of trap 
release in Eq. (12) [46]. This shows that all stationary quanti- 
ties of the released gas are described by a GGE. Furthermore, 
in Ref. [19] it has been shown that all non-equal time station- 
ary properties are always determined by the same ensemble 
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FIG. 4: The GGE structure factor S(k) as function of k/2kF 
{kp = 7rn) for different initial trap potentials ljN compared with 
the ground-state one (dashed line). 

describing the static quantities, and so, even in our case, they 
are encoded solely in the GGE. 

The structure factor in the GGE. The structure factor S{k) 
is the Fourier transform of the density-density correlation 
{n{x,t)n{0,t)) . In any ensemble which is diagonal in the 
Fourier modes, in the TD limit the structure factor can be writ- 
ten in terms of occupation modes as 

where the rhs is obtained using the GGE n^, given in Eq. (14). 
Here f[x) = [(4 + x^)E[l - A/x^) - ^{1 - 4/a;2)]|x|/6 
if |a:;| < 2 and zero otherwise where E{z) and K{z) are stan- 
dard elliptic functions and /(O) = 4/3. S{k) turns out to be 
an even function of k and monotonic for fc > 0. The plot 
of S{k) for different initial trapping potentials is reported in 
Fig. 4. S{k) resembles the one found numerically in [41] for 
the Lieb-Liniger gas. Because of the trap release constraint 
y/Nu > 2vln, we have S{k) > S{0) > 1 - S/Stt = 
0.151174. . . . This calculation shows how easy it is to ob- 
tain GGE predictions without solving the full non-equilibrium 
dynamics. 

The bosonic two-point function or one-body density matrix 
CB{x,y;t) = {^^x,t)^{y,t)) (with ^{y,t) bosonic anni- 
hilation operator) is a non-trivial quantity whose calculation 
presents difficulties also in thermal equilibrium [49]. How- 
ever, using the approach in [50], the computation is easy for 
large time and in the TD limit obtaining [46] 

Cs (x, y; i ^ oo) = C(x,y;i ^oo)e~2«|x-y|^ (jg^ 

with C{x, y;t —> oo) the fermion correlator in Eq. (12). For 
small distances, CB{x,y;t — > oo) shows a singular behavior 
of the form |a; — ?/| which is different from its thermal coun- 
terpart I a; — [49]. This behavior is strictly valid only in 
the TD limit because for any finite N, at very small distances 
Cb{x, y;t ^ oo) crosses over to \x — yp as expected from 
general arguments [49]. This finite N crossover is numeri- 
cally demonstrated in [46]. Consequently, the momentum dis- 
tribution function has a large momentum tail of the form fc~^ 



which crosses over to the standard k ^ for even larger k. This 
large-momentum crossover should be a measurable signature 
of the GGE. 

Trap to trap release. The case of a Tonks -Girardeau gas 
released not in a periodic system but in a larger harmonic trap 
has been solved by Minguzzi and Gangardt [36] who showed 
that the system oscillates forever without relaxation. How- 
ever, even in this case, it is simple to see that the time aver- 
aged two-point coiTelations (and hence by Wick theorem any 
observable) are still described by a GGE. 

Conclusions. In this letter we solved analytically the non 
equilibrium dynamics of a Tonks Girardeau gas following a 
trap release to a periodic geometry as in Fig. 1 . We prove that 
for long time and in the TD limit, any finite subsystem be- 
comes stationary and its behavior is described by a GGE. This 
provides the first analytic proof of a GGE for an inhomoge- 
neous initial state. We stress that the mechanism responsible 
for the equilibration is very different from the one in a global 
quantum quench since in the trap release it is due to the in- 
terference of the particles going around the circle many times. 
This equilibration mechanism is expected to be the same for 
any one dimensional gas released into a circle. 

Apart from the per se experimental interest [7, 8], these re- 
sults represent a first step towards a complete analytical un- 
derstanding of the famous quantum Newton cradle [2] at least 
in the Tonks Girardeau limit. 
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SUPPLEMENTARY MATERIAL 
The two-point fermionic correlation function 

Plugging the one-particle time-evolved wavefunctions [Eqs. (6) and (7) in the main text] in the definition of the fermionic 
correlator C{x, y;t) i^i t)^3 iVi 0' we have 

1 r uj^t 1 



p.q— — oo 

X 



/ x+pL \ f y + qL ^ 



Due to the oscillating phase factor, in the TD limit, the leading behavior of Eq. (S 1) is given only by the diagonal part p = q, so 
that 



C{x,y;t) ~ , > exp <^ i— — — — — f 1^ Xj { Xj , , (S2) 

that in the regime t ^ gives 

C,.„;,).l£e™£\,(i±f^),,(»±i^) 

p— — oo j— 

- . ATT. ° . /TT^ . /TT^ 



7V-1 



/oji -^-^ ^-^ \ \/ujt J \ \/u!t 

P— — 00 j—Q 



where we defined Xj{x) as the eigenfunctions of a harmonic oscillator with uj ~ 1 (i.e. Xj{x) — Xj(^)U=i)- 
In the large time limit t/L ^ oo (but with t/ L"^ ^ 0), the sum over p can be replaced by an integral 



L — 

where in the second line we used once again 7 — > 0. Now using 

^ N-l 

. [ 



3=0 



we have in the TD limit 



1 /-oc r°° rjh Ji V2ioN{x — y) 

C{x,y-t^oo) = - dfce'^(^-^)'=no(fc) = / — e^("-^)'=(vI/o|nfc|vI/o) = 2n 



i-Z-oo 27r ^/2u:N{x - y) 



(S4) 



(S5) 



where ^ — Lj {^/ojt) 0, allowing to rewrite C{x, y; t) as 

1 poo w-i 
C{x,y;t^^) = -j dze^v/^(--^)-^x,(7f +^)x,(7|+^) (S6) 

^ J —00 ■ r, 



(S7) 



lim I ^ |x,(a;)|' = %^2^^ ^ .^^^^^ ^^^^ 



(S9) 



which is Eq. (12) in the main text. It is also clear that this is nothing but the Fourier transform of (5'o|?^fc|^'o) obtained in Eq. 
(14) in the main text from the GGE, thus showing that infinite time limit and GGE results for C{x, y) are equal. 



(SIO) 



(Sll) 



The two-point bosonic correlation function 

The impenetrable bosons field operators $(x) are related to the fermionic ones ^{x) by the Jordan-Wigner transformation 

= exp|i7r J dz^\z)i'{z)^ ^{x), #^(a::) = <l^(a;) cxp |-i7r J (iz#^(z)#(z) 
Consequently, the equal time two-point bosonic correlation function Cb{x, y; t) with y > xis 

CB{x,y;t) = {^\x)^y)) = (^¥{x) cxp j-zTr^' dz*^(z)*(z)| ^(y)^ 
Expanding the exponential and using Wick theorem one has 

CB{x,y;t) = Y.^-^ / dzi--- dzr,{i>H^)^H^i)^{zi) ■ ■ ■ ^Hzn)-^{^n)^{y)) 

n=0 "'^ "'^ 

°° (—2)" 

= ^ j— / dzi--- dzn detC{xi,yj;t), 

n=0 "'^ "'^ 

where the indices i,j run from to n, and we fixed Xi — yi = Zi, Vi > 0, and xa = x, yo = y. The previous equation is a 
Fredholm's minor of the first order. 

In order to evaluate these correlators it is more convenient to introduce the N x N overlap matrix A{x, y; t) with elements 

rv 

A,jix, y;t)= dz'P*iz,t)'Pj{z,t), i, j & [0, . . . , N - 1]. (S12) 

J X 

in terms of which one has [50] 

N-l 

CB{x,y;t) = $*(a;,i)B,,(x,y;t)$,(2/,t), (S13) 

where the N x N matrix B(a;, y; t) is 

M{x,y;t) = det[¥]{p-^f, with P{x, y;t) = I - 2agn{y - x)A{x,y;t), (S14) 

and I is the N x N identity matrix. 

The large-time limit. Eq. (S 1 3) is a good starting point to evaluate analytically the large-time limit of the bosonic correlation 
function. Indeed, the one-particle time-evolved functions are given in Eqs. (6) and (7) of the main text, and proceeding as in Eq. 
(Si) we can write in the TD and large-time limits $*(2, t)<I>f,(z, t) as 

') - ^ E (^) « (^) - 11 .i^x: + .) » «) ^ i*..(si5, 

where in the last equality we used the orthonormality of the eigenfunctions Xa{x)- Consequently, the large-time behavior of the 
A and P matrices is 

Aab = J" dz^l{z,t)Mz,t) = ^^6at, P(a:,y;t^c»)= (^1-2^^^I. (S16) 

Clearly the above equations are valid as long as the rhs' are finite, i.e. when \x — y\/L ^ 0{1). For |a; — y| <C i different 
approaches must be used, as e.g. expanding the determinant in Eq. (Sll). 
From Eq. (S 16), the B matrix is 



l{x,y;t-^ oo) ^ fl - 2 ^"^ ^ I, and lim B(x, t -> cx)) = I e^^nla^-yl ^ (g^^-, 

y LI N~^oo 



where the large N limit has been taken keeping, as usual, n ~ N/L constant. Substituting Eq. (S17) in Eq. (S13), we finally 
have Eq. (16) in the main text, i.e. 

CB{x,y;t oo) = C(a;, y; t cx))e-2"l=^-''l. (S18) 
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FIG. 5: Exact bosonic correlation function Cb{x, y\t ^ oo) calculated by discretizing the Fredholm's minor in Eq. (SI 1). For large enough 
X, the data always agree with the prediction in Eq. (S18) (full lines), while for smaller x, the data approach it only for large enough uiN . The 
inset shows a zoom for very small x, for which the asymptotic \x\ behavior (dashed line) crosses over to a standard (non-singular) quadratic 
form. 



As anticipated in the main text, Eq. (S18) is valid only in the TD limit in the regime \x — y\/L ^ 0(1). For \x — y\ <^ L, 
the correlation function CB{x,y;t — > oo) crosses over to the standard singular behavior |a; — ?/|^, as expected from general 
arguments. In order to show the correctness of this statement, we calculate numerically Cb{x, 0; t — > oo) by discretizing the 
Fredholm's minor in Eq. (SI 1) as explained in Ref. [51] and using as input the GGE fermion correlation in Eq. (12) of the main 
text. In Fig. 5, we report the numerically calculated Cb{x, 0; t — oo) as function of x for different values of a; TV (we recall 
n = N/L is constant). It is clear that increasing N, the numerical data approach the asymptotic result in Eq. (S 18). However, if 
we zoom in the region of very small distances, as done in the inset of Fig. 5, the |a; — y\ singularity is absent, as expected, and 
the main singularity is of the form |a; — while the leading behavior is non-singular {x ~ y)'^. 



